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Abstract. - Logistic growth models are recurrent in biology, epidemiology, market models, and 
neural and social networks. They find important applications in many other fields including laser 
modelling. In numerous realistic cases the growth rate undergoes stochastic fluctuations and we 
consider a growth model with a stochastic growth rate modelled via an asymmetric Markovian 
dichotomic noise. We find an exact analytical solution for the probability distribution providing a 
powerful tool with applications ranging from biology to astrophysics and laser physics. 



Introduction. — In this letter we focus our attention on a growth model that was first 
used to describe the statistical behavior of a population of individual species; for example, 
human population growth. This is an area of interest several centuries old and is perhaps 
one of the oldest branches of biology to be studied quantitatively. The first model for human 
population growth was proposed by Malthus in 1798. In 1838 Vcrhulst 1 corrected this 
first model taking into account both the limitation of the growth of a population due to the 
competition between individuals, and the limitation on the density of the population that 
the environment can support. An example of this would be the limitation on the amount 
of food that the system is capable of producing. The proposed equation is known as the 
logistic equation, x = a x(l — x). 

There are many examples of assemblies that consist of a number of elements that interact 
through cooperative or competitive mechanisms. Some important examples include species 
that share a given environment, such as animals that live in the seas and rivers of our 
planet; the components of the central nervous system of a living being; transmission of 
diseases caused by different types of viruses; interacting vortices in a turbulent fluid; coupled 
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reactions between different chemical elements that make up our atmosphere; interactions 
between galaxies; competition between different political parties; business companies; and 
negotiating treaties of exchange between different countries. 

In the last few decades generalizations of the Malthus-Verhulst model have been applied 
to lasers physics [2"]I5] and have been widely considered in scientific literature. In Ref. [I], the 
reader can find a large list of references that use this type of model for many other physical 
processes such as the saturation growth process of a population [5HH] , which is considered one 
of the most successful models in the field of population dynamics. In addition the Malthus- 
Verhulst model has found applications the social sciences [9] [10], autocatalytic chemical 
reactions [IUH2], biological and biochemical systems [7], grain growth in polycrystalline 
materials [13] , cell growth in foam [14] , and as an effective model for the description of the 
populations of photons in a single mode laser [7] fT5Hl7] . The Lotka-Volterra model that was 
introduced early last century [IB] is a useful model to describe the interaction between two 
species. 

Mathematical models can be constructed either intuitively or from first principles to 
describe the phenomenon of competition and cooperation of many of the afore mentioned 
assemblies. This leads one to propose balance equations generally coupled and nonlinear. 
They contain some parameters that must be determined empirically or calculated from 
auxiliary equations. When the number of interacting variables is large, the number of 
balance equations is too large and therefore very difficult to solve. An example of this is 
classical mechanics applied to many-body systems. One does not generally know all the 
initial conditions and therefore it is necessary to develop statistical methods for multiple 
coupled rate equations that describe the behavior of the system far from equilibrium. Some 
important aspects of any assembly of elements that can be studied using statistical methods 
is its inherent stability, that is its stability with respect to small changes in growth rate and 
the introduction of new elements. 

The generalized Malthus-Verhulst model. — Our starting point is the generic 
stochastic differential equation for a growth model driven by a Markovian dichotomic noise 

x= [a (t)+a^{t)]x ^~ X ^ (1) 
M 

where the deterministic growth rate ao(t) is perturbed by the Markovian dichotomic noise 
in which a% and with fi > 0, are free parameters. Our main objective is to calculate 
the exact probability distribution for this generic model. 

The state space of the noise consists only of two levels, (Ai,— A2). This noise is 
called asymmetric Markovian dichotomic noise and it is also known as random telegraph 
noise. The temporal evolution of the conditional probability P(£,,t | £o>io) that completely 
characterizes the process is given by the following master equation [SIIT2] 



d_( P(A 1 ,t\^ th t ) \f-Ai A 2 \ / P(A u t\t ,t ) 
dt V -P(-A 2 ,t I 601*0) J \ Xi -A 2 J I P(-A 2 ,t I £o,*o) 



(2) 



where Ai and A2 are the probabilities by unit time of switching between states Ai and A 2 , 
so that Tj — 1/Xj arc the mean sojourn times in these states. The stationary solution of 
Eq. ([2]) can be obtained by setting 

P(£, 00 I 60, to) = - (X 2 Sa 1<6 + Ai<S_ A2 . 5 ) (3) 

7 

where 7 = Ai + A2. If the Markovian dichotomic noise has Eq. © as the initial condition, 
then £(t) is a stationary process. From Eq. ^ it follows that the mean value is 

<«*)> = A2Al " AlA2 - (4) 
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For the sake of simplicity, we require that the mean value (£(t)} vanish. This means that 

A 2 Ai = Ai A 2 = wo- (5) 

The correlation function is 

<£(*) £(*')> = ^ (Ai + A 2 ) 2 e-^-''l . (6) 
T 

Higher order correlation functions are more complicated. However, since the correlation 
function given by Eq. Q is indistinguishable from the Ornstein-Uhlenbeck process, the 
dichotomic noise found wide applications in building the models [5]. Furthermore, by an 
appropriate procedure of limit the dichotomic noise converges at the Gaussian white noise as 
the Ornstein-Uhlenbeck does, and it also converges at the white shot noise [15]. Experimen- 
tal evidences of the dichotomic noise have been found frequently in the litterature [5U2"0H2"2] . 

In Ref. [4] the author uses the direct method, which consists of formally integrating 
the stochastic differential equation and then taking the mean value over all realizations of 
the stochastic process. This method allows analytical treatment of the moments (x n (t)) 
for different types of noise, in particular for Gaussian white noise and white shot noise. 
In Ref. [23] the same author uses the inverse Mellin transform to calculate the stationary 
probability distribution. This last procedure appears to be difficult for two reasons. The 
first reason is that the mathematical problem of finding a distribution knowing its moments 
generally does not have a unique solution. The literature refers to this as the classical 
problem of moments [21] ■ The second reason is that indeed it is a very hard task to find 
an analytical inverse Mellin transform. For example if we consider the case with fj, = 1 in 
Eq. (P) and with £ (t) an Ornstein-Uhlenbeck process, the moments can be expressed as the 
following integral [25] 

(x n (t)) = ^=[ dz- - ^ ^n, n= 1,2,3,.... (7) 



in 



1 



(^) 



Vt) 



To find the inverse Mellin transform the parameter n has to be considered a real parameter. 
This fact makes it very difficult to perform the inversion. 

Using a relatively simple procedure the authors in Ref. [26 provide the exact probability 
distribution for a model like Eq. ([T]) when the noise £ (t) is given by the Ornstein-Uhlenbeck 
process. A further generalization can be found in Ref. [37] where several cases of white 
non-Gaussian noise are examined. 

For \i = 0, Eq. (TTJ reduces to the Gompertz model [5], x = [ao(t) + ai£_(t)] x\nx, and 
for fi = 1 it becomes the logistic equation driven by the dichotomic noise, namely 

i=(oo(t)+Bi£(t))x(l-ar). (8) 



Stochastic effects on Eq. ([8]) have been considered frequently in the literature. Refs. [5 H81fT8] 
contain several applications and developments of this model. In Refs. j28Tf3"Tj the transient 
behavior has been investigated when the system is driven by the same type of perturbation 
and the relaxation time of the system is calculated as a function of noise intensity. In 
Ref. [25] the results of Ref. [25J are extended to the case in which ao is perturbed by a 
colored Gaussian noise and confirmed by an analogical experiment, as well as by numerical 
simulations. To analyze a cancer cell population the authors of Ref. [32] consider the model 
x = ax—bx 2 +x^i(t)+^2(t) where and£ 2 (i) are correlated Gaussian white noises. They 
write the corresponding Fokker-Planck equation and analyze the behavior of the stationary 
probability density. 
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The exact probability distribution. — The model described by Eq. (JTJ) using the 
Stratonovich approach, can be reduced to an elementary differential equation by means of 
the transformation 



which leads to the equation 



" = '"1^-1 



y = -ao - ai£(i). (10) 



We emphasize that for the points x = and x = 1, the transformation Eq. (|9]) does not hold. 
The behavior of the system in these points has to be analyzed through a limit procedure. 
Following Rcf. 55], we can write the stochastic Liouville equation for the density function 
p(y, t; £) of a set of realizations of Eq. (ITU)) as 

d d c) 

—p(y,t;Z) = a —p(y,t;® + a 1 —Z(t)p(y,t-,Z). (11) 

Taking the mean value over all realizations of £ (t) we obtain 

dp dp dp 1 

d-t =a °d-y +ai ^y- (12) 

where p = p(y,t) = (p(y,t;Q) and p x = Pi(y,t) = (£(t)p(y,t;£)). Next, using the well- 
known Shapiro-Loginov formula for differentiation of exponentially correlated stochastic 
functions |34j . we obtain the following differential equation for the function p\(y,t) 

% =-A Pl + ao^+a 1 |(e 2 Wp( 2 /,^)) (13) 

where A = A x + A 2 . Following Ref. [33] we have £, 2 {t) = A 2 + A £(i), where A 2 = AiA 2 
and Aq = Ai — A 2 , transforms Eq. (fT3j) into 



&_-i, I + („ + «iV)^ + «A»* (14, 
Taking the time derivative of Eq. (TT2j) and combining it with Eq. (|14p we finally obtain 



d 2 p /q . a v d 2 P 
_-(2ao + a 1 A ) — 



<9 2 p 



dp .dp 



a oai A - a(A ) - a A^- + A-£ = 



ay 



'at 



(15) 



The Shapiro-Loginov formula has a hypothesis stating the statistical independence between 
£(t) and p(y,t;£) at the initial time t — 0. Therefore at the time i = 0, we have pi(y, 0) = 
(£(*)p(f)*)£)) |*=o= 0. Consequently the initial conditions for Eq. (fTS")) are 
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p(y,t)|t=o= 5 iv-va) 

The following change of variables 

t = T, y = z-a T 
further simplifies Eq. (|15j) . so that we end up with 

d 2 



p{y,t)\t=o= a S' (y - y ) ■ 



— -a A 92 
dr 2 1 ° drdz 



dz 



dr 



P(z,t) = 



(16) 
(17) 

(18) 



where P(z,t) = p(y,t). A formal solution of Eq. (|T5)) . satisfying the initial conditions (TTJ 
and vanishing for z — > ±oo, is given by 
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P(z,t) = e"* T 



ikr 



A A — iai Aofc . A 

cos — r H sm — r 

2 A 2 



2^ 



(19) 



where 



aiA 



-T + Z — Zq 



Evaluating the integral we obtain 



, A = Y (Ai + A 2 ) 2 a\k 2 + 2ia 1 A Ak - A 2 . 



(20) 



P(z,t) 



A — aiAoM 



-9(vt- |r|) 



ttT + Ol A £ h 

. J 1 

yj V 2 T 2 — T 2 



la (by/l 



i . aiAo \ . , , ( QiA 
1 H — 1 o (vt — r) + 1 — I o (vt + r) 



2v 



2v 



(21) 



where r, as function of z, is defined by Eq. (I2TJ1) . /„((?) are the modified Bessel's functions, 
9{q) is the step function, 5{q) is the Dirac's delta. Furthermore, we defined in a compact 

manner 

aiAA , % /ATA^A|ai| 



b = 



Ai + A 2 



With the help of Eqs. 
variables a?, i as 



4v 2 ' 2w 2 2 

|) and (|T7)) we may write the solution as function of the original 



p(x,t) 



x(l- x^) 



P 



In 



1-x" 

jJLX^ 



agt, t 



(22) 



Numerical Analysis. — In the following we numerically implement the dynamic equa- 
tion (ITUl) driving the diffusion of the variable y. An ensemble of N "trajectories" of the vari- 
able £(t), the dichotomic noise in (fTOJ) . is produced via a generator of random numbers with 
poisson distribution which returns the time intervals in which the stochastic dichotomous 
variable £(t) retains either of its two values. For each trajectory of £(t), a trajectory of the 
variable y is obtained, integrating Eq. (fTU|) . The probability for a given value of y at time 
t is then calculated as a simple average over the ensemble of N trajectories so obtained. A 
subsequent conversion to the x-space through the transformation ^ allows us to obtain the 
probability density P(x,t). The dichotomous process £(t) is assumed in a stationary condi- 
tion at time t = 0, i.e. at such a time A^Ai/(Ai + A 2 ) trajectories are taken with initial value 
Aj for the variable £(t) and NX 2 /(Xi + A 2 ) with the value — A 2 . The numerical results are 
compared to the analytical expression given by Eqs. (|2"Tj) and (|2"2")l in Fig. (fT]) which shows 
a perfect agreement. 

Analysis of the results. — Thanks to calculations performed in the previous sections 
we have at our disposal the exact solution for any value of the parameter fi of physical 
interest. A detailed study of P(x, t) is beyond the purpose of this letter. We shall limit 
ourselves to rediscover some results known in the literature for the case \i = 1, which is the 
Malthus-Verhulst model, and for fj, = 0, which is the Gompertz model. 

We note a first result of solution (fTSj) . specifically of Eq. (|2Tj) . Due to the asymmetric case 
under consideration, there exists a particular choice of parameters such that a coefficient of 
the two deltas vanishes. To fix the ideas let us set a\ > 0, then we can make one of the 
two delta coefficients vanish if aiAo = 2v, that is to say A 2 = 0. Note that A 2 = does 
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Fig. 1: The probability distribution at times t = 5.7 (right curve) and t = 7.7 (left curve) with 
parameters Ai = 1.2, A2 = 2.8, ao = 1.0, a\ = 1.0, fi — 1.0, Ai = 0.6, A2 = 1.4 . The diamonds are 
the result of the numerical simulation, the continuous curves are the analytic solutions as given by 
Eqs. (|2ip and (|22|l . The dirac deltas have been subtracted in both cases. 



not imply a vanishing value of the parameter v that represents the propagation speed of the 
peaks. 

Taking appropriate limits on the parameters A, Ai, and A2 we can rediscover different 
well known stochastic processes. Following Ref. [19], we consider the limit 

Ai = A2 = A — > 00, Ai = A2 = A — > 00 

which corresponds to the Gaussian white noise. Keeping constant the ratio A 2 /A we obtain 



P(z,r) 



,-At 



-At 



1 (vt— I r\) 



bl\ ( bvr 

1 2vt 



-To ( i"'T 



br 2 
2vt 



2v 



5 (vt — r) + 1 



2v 



5 (vt 



(23) 



Using the asymptotic expression for the Bessel's functions that is I n (x) ~ exp[x]/v27ra' for 
x — > 00, we finally obtain 



P(z,t) 



1 



V2vrL>r 



exp 



2Dt 



(24) 



where by definition D — a 2 A 2 / A and we neglect the two exponentially damped deltas. 

Still following Ref. [H] , to obtain the white shot noise limit we have to take the symmetric 
dichotomous noise limit 



Ai = A 2 = A -)• 00, Ai = A 2 = A 



A 

00, — 



7 



where the 7 parameter is called non-Gaussianity parameter. The relation with the diffusion 
coefficient D is given by D = j 2 \. In the above limit bv = A — > 00 so that we again end up 
with Eq. (|23| because the argument of Bessel's functions becomes infinite. Consequently we 
rediscover Eq. (|2"4"|) 

As the last result we consider the limit for /x — > which leads us to the Gompertz model. 
The transformation (El) becomes 
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y = ln[— In x\. 



(25) 



For brevity we consider the symmetric case so that A = 0. Using Eq. (|2"2"j) we find that the 
asymptotic solution is 



P(s,t)«- 



1 



1 



y/2nDt x\nx 



exp 



(ln[— In a;] + ao<) 2 
2Dt 



(26) 



An analysis of this solution shows that P(x,i), given by Eq. ([26| . always has a minimum 
located near the origin and a maximum located near x — 1. Finally P(x, t) diverges at x = 
as 

P(x,t) ' ' 



(— In a;) 



These results are graphed in Fig. [2] 



P(x,l) 




X 



Fig. 2: The probability distribution P(x,t) for the Gompertz model at times t = 1. The value of 
the parameters are ao — 1 and D = 1. 



Concluding remarks. — In this letter we considered a growth model with a stochastic 
growth rate modelled with a Markovian dichotomic noise. We found an exact solution for 
the generic parameters /z, Ai and A2 which makes it possible to apply the result to a variety 
of models ranging from biology [5J[TH] to astrophysics [35] and laser physics. For example 
our model fits the the dichotomic noise condition found in ion channel experiments |38f - 
HD] . Patch-clamp experiments detect a fluctuating dichotomous current which is related to 
the number of ions flowing in the channel. Yet, in single mode lasers, the growth model 
with [i = 2 and a constant growth rate, describes the time evolution for the N th mode 
of the electric field in the laser cavity (2][3]. The model in this letter makes a case for 
dichotomic fluctuations around a resonance frequency, providing opportunities for additional 
applications. An interpretation of the Malthus-Verhulst model which describes the photon 
population in a laser cavity, corresponding to the case \i = 1, was introduced by McNeil 
and Walls [37]. We found these to be interesting cases and chose to examine the Malthus- 
Verhulst model, which has a very wide range of applications such as dynamics of populations, 
chemical reactions, and photon population in a laser cavity [TT | I37 ] . I n addition we studied 
the Gompertz model, corresponding to /1 = 0, when slightly modified, is used to describe 
the tumor growth dynamics. We confirmed our results with a numerical simulation showing 
perfect agreement with the analytical formulas. In this letter we focused only on two models, 
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but the obtained results may applied to other systems that will be the subject of future 
research. 
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